%     VIRTUAL FABRIC - программа для получения эффективных свойств текстильных композитов
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function [NODES, ELEMENTS, ORIENTATIONS, RVE_Coordinates, Node_Merge_Tolerance, Fout, V_resin_free] = Process_Voxel_Model(Voxel_Model_Directory, Voxel_Model_Name, Work_Directory, Job_Name, Analysis_Type, Part_Name, Mechanical_Element_Type)

global ErrCode;
global Merge_Tolerance_Multiplier;

NODES = [];
ELEMENTS = [];
ORIENTATIONS = [];
RVE_Coordinates = [];
Node_Merge_Tolerance = [];
Fout = [];
V_resin_free = [];

% Имя модели
Model_Name = Job_Name;

% Путь к файлу с воксельной моделью
File_path_model = [Voxel_Model_Directory '\' Voxel_Model_Name '.inp'];

% Открытие файла с воксельной моделью
Fin_model = fopen(File_path_model, 'r');

if (Fin_model < 0)
    fprintf (1, '\nОшибка функции Process_Voxel_Model: не удалось открыть файл [%s] !\n', File_path_model);
    ErrCode = -1;
    return;
end

% Путь к файлу с ориентацией
File_path_ori = [Voxel_Model_Directory '\' Voxel_Model_Name '.ori'];

% Открытие файла с ориентацией
Fin_ori = fopen(File_path_ori, 'r');

if (Fin_ori < 0)
    fprintf (1, '\nОшибка функции Process_Voxel_Model: не удалось открыть файл [%s] !\n', File_path_ori);
    ErrCode = -1;
    fclose(Fin_model);
    return;
end

% Путь к выходному файлу
File_path_out = [Work_Directory '\' Job_Name '.inp'];

% Открытие выходного файла
Fout = fopen(File_path_out, 'wt');

if (Fout < 0)
    fprintf (1, '\nОшибка функции Process_Voxel_Model: не удалось открыть для записи файл [%s] !\n', File_path_out);
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    return;
end

% Начало выходного файла
fprintf (Fout, '*Heading\n');

fprintf (Fout, '** Job name: %s Model name: %s\n', Job_Name, Model_Name);
fprintf (Fout, '** Generated by: Virtual_Fabric\n');

fprintf (Fout, '*Preprint, echo=NO, model=NO, history=NO, contact=NO\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** PARTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Part, name=%s\n', Part_Name);

% Чтение файла с воксельной моделью

Nodes_alloc_param = 500000;
Elm_alloc_param = 500000;
Ori_alloc_param = 500000;

Node_counter = 0;
Elm_counter = 0;
Ori_counter = 0;

NODES = zeros(Nodes_alloc_param, 4);
ELEMENTS = zeros(Elm_alloc_param, 9);
ORIENTATIONS = zeros(Ori_alloc_param, 7);

Node_Stop_Flag = 0;
Set_Stop_Flag = 0;

Matrix_Elm_Counter = 0;

DISTRIBUTION_TABLE_Flag = 0;
DISTRIBUTION_Flag = 0;
ORIENTATION_Flag = 0;
SOLID_SECTION_Flag = 0;

while (~feof(Fin_model))
    
    string = strtrim(fgets(Fin_model));
    
    if (isempty(string))
        continue;
    end
    
    if (string(1) ~= '*')
        continue;
    end
    
    if ((length(string) >= 2) && (string(1) == '*') && (string(2) == '*'))
        
        if (~isempty(strfind(string, 'PERIODIC BOUNDARY CONDITIONS'))) %#ok<STREMP>
            Node_Stop_Flag = 1;
            Set_Stop_Flag = 1;
        end
        
        continue;
    end
    
    temp_cell_array = strsplit(string(2:length(string)), {','});
    
    if (isempty(temp_cell_array{1}))
        fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
        ErrCode = -1;
        fclose(Fin_model);
        fclose(Fin_ori);
        fclose(Fout);
        Fout = [];
        return;
    end
    
    switch upper(strtrim(temp_cell_array{1}))
        
        case 'NODE'
            
            if (Node_Stop_Flag == 1)
                continue;
            end
            
            fprintf (Fout, '%s\n', string);
            
            while(~feof(Fin_model))
                
                string = fgets(Fin_model);
                strlen = length(string);
                string = strtrim(string);
                
                if (isempty(string))
                    break;
                end
                
                if (string(1) == '*')
                    fseek(Fin_model, -strlen, 'cof');
                    break;
                end
                
                [tmp, count] = sscanf(string, '%e,');

                if (count == 4)
                    
                    Node_counter = Node_counter + 1;
                    
                    if (Node_counter > size(NODES,1))
                        NODES(size(NODES,1) + Nodes_alloc_param,1) = 0;
                    end
                    
                    NODES(Node_counter,:) = tmp';
                    
                    if (Node_counter == 1)
                        
                        RVE_Coordinates = zeros(3,2);
                        RVE_Coordinates(1,1) = NODES(Node_counter,2);
                        RVE_Coordinates(1,2) = NODES(Node_counter,2);
                        RVE_Coordinates(2,1) = NODES(Node_counter,3);
                        RVE_Coordinates(2,2) = NODES(Node_counter,3);
                        RVE_Coordinates(3,1) = NODES(Node_counter,4);
                        RVE_Coordinates(3,2) = NODES(Node_counter,4);
                        
                    else
                        
                        if (NODES(Node_counter,2) < RVE_Coordinates(1,1))
                            RVE_Coordinates(1,1) = NODES(Node_counter,2);
                        end
                        if (NODES(Node_counter,2) > RVE_Coordinates(1,2))
                            RVE_Coordinates(1,2) = NODES(Node_counter,2);
                        end
                        
                        if (NODES(Node_counter,3) < RVE_Coordinates(2,1))
                            RVE_Coordinates(2,1) = NODES(Node_counter,3);
                        end
                        if (NODES(Node_counter,3) > RVE_Coordinates(2,2))
                            RVE_Coordinates(2,2) = NODES(Node_counter,3);
                        end
                        
                        if (NODES(Node_counter,4) < RVE_Coordinates(3,1))
                            RVE_Coordinates(3,1) = NODES(Node_counter,4);
                        end
                        if (NODES(Node_counter,4) > RVE_Coordinates(3,2))
                            RVE_Coordinates(3,2) = NODES(Node_counter,4);
                        end
                        
                    end
                    
                else
                    fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
                    ErrCode = -1;
                    fclose(Fin_model);
                    fclose(Fin_ori);
                    fclose(Fout);
                    Fout = [];
                    return;
                end
                
                fprintf (Fout, '%s\n', string);
                
            end
            
        case 'ELEMENT'
            
            if (numel(temp_cell_array) < 2)
                fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
                ErrCode = -1;
                fclose(Fin_model);
                fclose(Fin_ori);
                fclose(Fout);
                Fout = [];
                return;
            end
            
            if (isempty(temp_cell_array{2}))
                fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
                ErrCode = -1;
                fclose(Fin_model);
                fclose(Fin_ori);
                fclose(Fout);
                Fout = [];
                return;
            end
            
            switch upper(strtrim(temp_cell_array{2}))
                
                case 'TYPE=C3D8'
                    % Все нормально
                    
                case 'TYPE=C3D8R'
                    % Все нормально
                    
                otherwise
                    
                    fprintf (1, '\nОшибка: недопустимый тип элементов в файле с воксельной моделью в строке [%s] !\n', string);
                    ErrCode = -1;
                    fclose(Fin_model);
                    fclose(Fin_ori);
                    fclose(Fout);
                    Fout = [];
                    return;
                
            end
            
            switch Analysis_Type
                
                case {'MECHANICAL', 'MESH_GENERATOR'}
                    
                    switch Mechanical_Element_Type
                        case 'C3D8R'
                            fprintf (Fout, '*Element, type=C3D8R\n');
                        case 'C3D8RH'
                            fprintf (Fout, '*Element, type=C3D8RH\n');
                        case 'C3D8'
                            fprintf (Fout, '*Element, type=C3D8\n');
                        case 'C3D8H'
                            fprintf (Fout, '*Element, type=C3D8H\n');
                        case 'C3D8I'
                            fprintf (Fout, '*Element, type=C3D8I\n');
                        case 'C3D8IH'
                            fprintf (Fout, '*Element, type=C3D8IH\n');
                    end
                    
                case 'HEAT_TRANSFER'
                    fprintf (Fout, '*Element, type=DC3D8\n');
                
            end
            
            while(~feof(Fin_model))
                
                string = fgets(Fin_model);
                strlen = length(string);
                string = strtrim(string);
                
                if (isempty(string))
                    break;
                end
                
                if (string(1) == '*')
                    fseek(Fin_model, -strlen, 'cof');
                    break;
                end
                
                [tmp, count] = sscanf(string, '%i,');

                if (count == 9)
                    
                    Elm_counter = Elm_counter + 1;
                    
                    if (Elm_counter > size(ELEMENTS,1))
                        ELEMENTS(size(ELEMENTS,1) + Elm_alloc_param,1) = 0;
                    end
                    
                    ELEMENTS(Elm_counter,:) = tmp';
                    
                else
                    fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
                    ErrCode = -1;
                    fclose(Fin_model);
                    fclose(Fin_ori);
                    fclose(Fout);
                    Fout = [];
                    return;
                end
                
                fprintf (Fout, '%s\n', string);
                
            end
            
        case 'DISTRIBUTION TABLE'
            
            DISTRIBUTION_TABLE_Flag = 1;
            
            string = strtrim(fgets(Fin_model));
            
            temp_cell_array = strsplit(string(1:length(string)), {','});
            
            if (numel(temp_cell_array) ~= 2)
                fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
                ErrCode = -1;
                fclose(Fin_model);
                fclose(Fin_ori);
                fclose(Fout);
                Fout = [];
                return;
            end
            
            if (isempty(temp_cell_array{1}) || isempty(temp_cell_array{2}))
                fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
                ErrCode = -1;
                fclose(Fin_model);
                fclose(Fin_ori);
                fclose(Fout);
                Fout = [];
                return;
            end
            
            if (~strcmpi(strtrim(temp_cell_array{1}), 'COORD3D'))
                fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
                ErrCode = -1;
                fclose(Fin_model);
                fclose(Fin_ori);
                fclose(Fout);
                Fout = [];
                return;
            end
            
            if (~strcmpi(strtrim(temp_cell_array{2}), 'COORD3D'))
                fprintf (1, '\nОшибка: нарушен формат файла c воксельной моделью в строке [%s] !\n', string);
                ErrCode = -1;
                fclose(Fin_model);
                fclose(Fin_ori);
                fclose(Fout);
                Fout = [];
                return;
            end
            
        case 'DISTRIBUTION'
            
            DISTRIBUTION_Flag = 1;
            
            fprintf (Fout, '*Distribution, name=Orientation_Vectors, location=ELEMENT, Table=Orientation_Vectors_Table\n');
            
            First_Ori_String_Flag = 1;
            
            while(~feof(Fin_ori))
                
                string = strtrim(fgets(Fin_ori));

                if (isempty(string))
                    continue;
                end
                
                if (string(1) == '*')
                    continue;
                end
                
                if (First_Ori_String_Flag == 0)
                    
                    [tmp, count] = sscanf(string, '%e,');

                    if (count == 7)

                        Ori_counter = Ori_counter + 1;

                        if (Ori_counter > size(ORIENTATIONS,1))
                            ORIENTATIONS(size(ORIENTATIONS,1) + Ori_alloc_param,1) = 0;
                        end

                        ORIENTATIONS(Ori_counter,:) = tmp';

                    else
                        fprintf (1, '\nОшибка: нарушен формат файла c ориентацией в строке [%s] !\n', string);
                        ErrCode = -1;
                        fclose(Fin_model);
                        fclose(Fin_ori);
                        fclose(Fout);
                        Fout = [];
                        return;
                    end
                    
                end
                
                fprintf (Fout, '%s\n', string);
                
                if (First_Ori_String_Flag == 1)
                    First_Ori_String_Flag = 0;
                end
                
            end
            
        case 'ORIENTATION'
            
            ORIENTATION_Flag = 1;
            
            fprintf (Fout, '*Orientation, name=Ori_Orientation_Vectors, DEFINITION=COORDINATES, system=RECTANGULAR\n');
            fprintf (Fout, 'Orientation_Vectors\n');
            fprintf (Fout, '1, 0.\n');
            
        case 'ASSEMBLY'
            
            Node_Stop_Flag = 1;
            Set_Stop_Flag = 1;
            
        case 'ELSET'
            
            if (Set_Stop_Flag == 1)
                continue;
            end
            
            Matrix_Set_Flag = 0;
            
            if ((numel(temp_cell_array) >= 2) && strcmpi(strtrim(temp_cell_array{2}), 'ElSet=Matrix'))
                Matrix_Set_Flag = 1;
            end
            
            string = strrep(string, 'AllElements', 'ALL_ELEMENTS');
            
            fprintf (Fout, '%s\n', string);
            
            while(~feof(Fin_model))
                
                string = fgets(Fin_model);
                strlen = length(string);
                string = strtrim(string);
                
                if (isempty(string))
                    break;
                end
                
                if (string(1) == '*')
                    fseek(Fin_model, -strlen, 'cof');
                    break;
                end
                
                if (Matrix_Set_Flag == 1)
                    [tmp, count] = sscanf(string, '%i,'); %#ok<ASGLU>
                    Matrix_Elm_Counter = Matrix_Elm_Counter + count;
                end
                
                fprintf (Fout, '%s\n', string);
                
            end
            
        case 'NSET'
            
            if (Set_Stop_Flag == 1)
                continue;
            end
            
            string = strrep(string, 'AllNodes', 'ALL_NODES');
            
            fprintf (Fout, '%s\n', string);
            
            while(~feof(Fin_model))
                
                string = fgets(Fin_model);
                strlen = length(string);
                string = strtrim(string);
                
                if (isempty(string))
                    break;
                end
                
                if (string(1) == '*')
                    fseek(Fin_model, -strlen, 'cof');
                    break;
                end
                
                fprintf (Fout, '%s\n', string);
                
            end
            
        case 'SOLID SECTION'
            
            SOLID_SECTION_Flag = 1;
            
            string = strrep(string, 'TexGenOrientations', 'Ori_Orientation_Vectors');
            string = strrep(string, 'Mat0', 'Resin_Material');
            string = strrep(string, 'Mat1', 'Yarn_Material');
            
            fprintf (Fout, '%s\n', string);
            
            while(~feof(Fin_model))
                
                string = fgets(Fin_model);
                strlen = length(string);
                string = strtrim(string);
                
                if (isempty(string))
                    break;
                end
                
                if (string(1) == '*')
                    fseek(Fin_model, -strlen, 'cof');
                    break;
                end
                
                fprintf (Fout, '%s\n', string);
                
            end
        
    end
    
end

fprintf (Fout, '*End Part\n');

fprintf (Fout, '*Distribution Table, Name=Orientation_Vectors_Table\n');
fprintf (Fout, 'coord3d, coord3d\n');

% Проверка флагов
if (DISTRIBUTION_TABLE_Flag == 0)
    fprintf (1, '\nОшибка: в файле с воксельной моделью нет записи *DISTRIBUTION TABLE !\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
end

if (DISTRIBUTION_Flag == 0)
    fprintf (1, '\nОшибка: в файле с воксельной моделью нет записи *DISTRIBUTION !\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
end

if (ORIENTATION_Flag == 0)
    fprintf (1, '\nОшибка: в файле с воксельной моделью нет записи *ORIENTATION !\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
end

if (SOLID_SECTION_Flag == 0)
    fprintf (1, '\nОшибка: в файле с воксельной моделью нет записей *SOLID SECTION !\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
end

% Финализация массивов
if (Node_counter == 0)
    fprintf (1, '\nОшибка: в файле с воксельной моделью нет узлов !\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
else
    NODES = NODES(1:Node_counter,:);
end

if (Elm_counter == 0)
    fprintf (1, '\nОшибка: в файле с воксельной моделью нет элементов !\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
else
    ELEMENTS = ELEMENTS(1:Elm_counter,:);
end

if (Ori_counter == 0)
    fprintf (1, '\nОшибка: в файле с ориентацией нет векторов ориентации !\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
else
    ORIENTATIONS = ORIENTATIONS(1:Ori_counter,:);
end

% Проверка нумерации узлов
for i=1:size(NODES,1)
    if (NODES(i,1) ~= i)
        fprintf (1, '\nОшибка функции Process_Voxel_Model: некорректная нумерация узлов!\n');
        ErrCode = -1;
        fclose(Fin_model);
        fclose(Fin_ori);
        fclose(Fout);
        Fout = [];
        return;
    end
end

% Проверка нумерации элементов
for i=1:size(ELEMENTS,1)
    if (ELEMENTS(i,1) ~= i)
        fprintf (1, '\nОшибка функции Process_Voxel_Model: некорректная нумерация элементов!\n');
        ErrCode = -1;
        fclose(Fin_model);
        fclose(Fin_ori);
        fclose(Fout);
        Fout = [];
        return;
    end
end

% Проверка нумерации векторов ориентации
for i=1:size(ORIENTATIONS,1)
    if (ORIENTATIONS(i,1) ~= i)
        fprintf (1, '\nОшибка функции Process_Voxel_Model: некорректная нумерация векторов ориентации!\n');
        ErrCode = -1;
        fclose(Fin_model);
        fclose(Fin_ori);
        fclose(Fout);
        Fout = [];
        return;
    end
end

% Проверка совпадения размеров массива элементов и массива векторов ориентации
if (size(ORIENTATIONS,1) ~= size(ELEMENTS,1))
    fprintf (1, '\nОшибка функции Process_Voxel_Model: размеры массива элементов и массива векторов ориентации не совпадают!\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
end

% Нахождение длины минимального ребра элемента
Min_Element_Edge = min([norm(NODES(ELEMENTS(1,2), 2:4) - NODES(ELEMENTS(1,3), 2:4)), ...
                        norm(NODES(ELEMENTS(1,2), 2:4) - NODES(ELEMENTS(1,5), 2:4)), ...
                        norm(NODES(ELEMENTS(1,2), 2:4) - NODES(ELEMENTS(1,6), 2:4))]);

% Допуск для слияния узлов
Node_Merge_Tolerance = Min_Element_Edge * Merge_Tolerance_Multiplier;

% Проверка числа элементов матрицы
if (Matrix_Elm_Counter == 0)
    fprintf (1, '\nОшибка функции Process_Voxel_Model: в воксельной модели отсутствует матрица!\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
end

if (Matrix_Elm_Counter == size(ELEMENTS,1))
    fprintf (1, '\nОшибка функции Process_Voxel_Model: в воксельной модели отсутствуют жгуты!\n');
    ErrCode = -1;
    fclose(Fin_model);
    fclose(Fin_ori);
    fclose(Fout);
    Fout = [];
    return;
end

% Объемная доля свободного связующего
V_resin_free = Matrix_Elm_Counter / size(ELEMENTS,1);

% Закрытие входных файлов
fclose(Fin_model);
fclose(Fin_ori);
